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^ , Abstract. We investigate the system size scaling of the net defect number 

a^^ ■ created by a rapid quench in a second-order quantum phase transition from 

' an 0{N) symmetric state to a phase of broken symmetry. Using a controlled 

mean-field expansion for large TV, we find that the net defect number variance 
in convex volumina scales like the surface area of the sample for short- 
range correlations. This behaviour follows generally from spatial and internal 
C^ . symmetries. Conversely, if spatial isotropy is broken, e.g., by a lattice, and in 

addition long-range periodic correlations develop in the broken-symmetry phase, 
C^ ' we get the rather counterintuitive result that the scaling strongly depends on 

the dimension being even or odd: For even dimensions, the net defect number 
variance scales like the surface area squared, with a prefactor oscillating with the 
i-pj , system size, while for odd dimensions, it essentially vanishes. 
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1. Introduction 

Rapidly sweeping through a second-order (quantum or thermal) symmetry-breaking 
phase transition in general may trigger the creation of topological defects. This is due 
to the fact that the system cannot adjust abruptly to the new order, so that (time- 
dependent) correlations are spontaneously generated in the interacting many-body 
system. Correctly counting the topological defects generated by the sweep amounts 
to an analysis of these correlations. 

The fundamental global requirement for topological defect generation is that the 
symmetry group of the new phase permits the defects in the sense that the homotopy 
group of the coset space, formed by the quotient space of unbroken (old) and broken 
(new) symmetries is nontrivial. The physical consequences of the latter mathematical 
fact were first used by Kibble [1] in the high-energy context of symmetry breakings 
from a (grand unified) symmetry group of generally high dimension to one of lower 
dimension in the earliest stages of the universe. This led to the classification of the 
symmetry group (s) for which cosmic strings and domain walls can, in principle, be 
produced by the universe's expansion after the big bang. An implementation of this 
idea directly accessible to laboratory-scale experimental realizations in condensed- 
matter systems was devised nine years later by Zurek |2]. Hence the historically 
developed and by now established name Kibble-Zurek mechanism for this particular 
symmetry-breaking based scenario of topological defect formation. While Kibble 
realized that in the early universe relativistic causality alone mandates the appearance 
of defects according to the nontrivial homotopy groups relevant to the symmetry- 
breaking scheme at hand, Zurek used scalings of the relaxation time and healing 
length with the quench time in the vicinity of the critical point to predict topological 
defect densities. The scaling behaviors derived by Zurek, based on the critical 
exponents of the phase transition under consideration, therefore led to concretely 
testable experimental signatures of the general idea of Kibble. 

The basic physical phenomenon is depicted in Fig.[T]for the simple case of a U{\) 
order parameter in the plane: It relies on the finite propagation speed communicating 
the new order and the corresponding correlations between spatially well-separated 
regions. These finite communication speeds lead to the whole sample not being able 
to establish the new order homogeneously throughout, thus causing an inhomogeneous 
order parameter distribution. Order parameter domains form, which, with some 
appropriate statistical probability, potentially enclose topological defects. 

The Kibble-Zurek mechanism may occur in many diverse physical settings, for 
instance, in nonequilibrium phase transitions during the early universe [31 E], and 
has been measured in various condensed matter systems like liquid crystals [5] or 
superconducting films |6]. There has been a tremendously increasing interest in the 
closely related nonequilibrium quantum many-body phenomena during quantum phase 
transitions in recent years, cf., e.g., [7l|8l[9llIll|ll[12l[13l[14l[15l[16l[17l[18l[TS|23 
mj [55] , due also to the fact that for controlled dynamical entanglement generation 
[23j the dynamical evolution of quantum many-body correlations is strongly relevant. 

A clean experimental evidence for Kibble-Zurek mediated defect formation in 
dilute cold atomic gases, which are an experimentally precisely in situ controllable 
variant of conventional condensed matter systems, has been obtained in spinor as 
well as scalar Bose gases [241 [55] . These dilute gases are particularly suitable to 
test predictions of the dynamical Kibble-Zurek theory because of their comparatively 
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long re-equilibration times, which effectively permit a real-time study of the evolution 
of perturbations in the spatial many-body correlations and thus ultimately also a 
study of the creation process of topological defects. Furthermore, it was shown that 
the Kibble-Zurek type scaling of the defect number with the sweep rate through the 
phase transition depends on the inhomogenity of the sample [^ \T7\ (is becoming 
more strongly dependent on the sweep rate), which is relevant for the intrinsically 
inhomogenous harmonically trapped gases. While most investigations to date have 
been carried out for bosonic systems in low dimension (where dimension here refers 
to real as well as order-parameter dimension), due to the rapidly growing interest in 
cold atomic gases, more recently fermionic problems in these systems have also been 
tackled, e.g. in [281129] . 

When calculating defect densities from correlation functions, frequently the 
excitation densities (quasiparticle densities) and the topological defect densities 
actually created by the quench are assumed to be in a one-to-one correspondence 
cf., e.g., [IHlIiniHT]. This identification is however misleading in general, because 
the necessary correlations created by the quench are much stronger for creating 
topological defects as compared to the correlations for quasiparticle excitations in 
high dimension [30]. For an illustration of this fact cf. Figl2l A topological defect at a 
given vertex enclosed by domains after the quench implies that the direction vectors 
of all neighbouring domains either point away or towards the vertex (for a monopole 
configuration of the spin vectors, see also below). In a high-dimensional space, for 
example, this condition for a defect core to be enclosed by the domains at the vertex 
such that a winding of the order parameter orientation is induced becomes increasingly 
difficult to fulfill, given that the direction vectors orient themselves statistically inside 
the sample (i.e. due to quantum or thermal fluctuations establishing the order 




Figure 1. Origin of topological defect formation due to the Kibble- 
Zurek mechanism: Rapidly sweeping through a second-order phase transition 
statistically generates domains due to the finite propagation speed communicating 
of the new order in space and establishing the corresponding order parameter 
correlations. If by chance, at a given vertex, the domains meet such that a 
(generalized) "phase winding" is produced, a topological defect with a core of 
the old symmetry (yellow) is formed. In the simple example shown in the figure, 
a {/(I) order parameter in the two-dimensional plane with field direction indicated 
by the red arrows, this amounts to the order parameter vector orientation winding 
around the core by 27r, for the planar spin hedgehog depicted. 
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Figure 2. In high dimensions, a rapidly increasing degree of coordination of 
the order parameter orientation in neighboring domains is necessary to enclose a 
defect. In the three-dimensional example shown above for the classic monopole 
configuration, all arrows in the surrounding domains (red) have to point away 
from the central defect core (yellow), cf. the two-dimensional situation in Fig.[T] 
which already requires significantly less coordination. 



parameter correlations) . 

In the following, we will investigate topological defect creation in quantum 
quenches concentrating on the behaviour of the defect creation rate for high internal 
and spatial dimension of the system. To this end, we derive essentially exact 
statements on topological defect production in the large N limit of an 0{N) symmetric 
model in _D = iV spatial dimensions. Similar to large- A'' approaches in condensed 
matter and field theory we argue that our predictions for the defect production will 
remain qualitatively valid also for smaller dimensions. More specifically, assuming that 
there is no critical value of N where the system behaviour changes in a discontinous 
manner, we expect our results to be qualitatively applicable also to finite N, e.g., 
A'^ = 3, which are accessible to experimental tests, for example in spinor Bose-Einstein 
condensates with large spin, cf. the review |31| . 

In the next Section, we begin by deriving a large N expansion for a general 0{N) 
model resulting in a bilinear effective action, which is used to obtain the equations of 
motion describing the instability of direction modes during the phase transition. We 
then apply these equations to calculate the behavior of the spatial direction correlation 
functions shortly after the transition. The latter result is used to deduce ab initio the 
winding number variance for our 0{N) model. Finally, we delineate the behavior 
of the winding number variance with the size of the sample in dependence on the 
spatial behavior of the direction correlations (which is, e.g., short-range or long-range 
periodic). 

2. Large- A^ Expansion 

In order to illustrate the main idea - i.e., how to obtain an effective quadratic action 
for a iV-component field $ = ($i, $2, ■ • ■ 7 ^n) in the large N limit, let us consider 
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the (Euclidean) generating functional of the non- linear 0{N) model 

Zp[J]= f ^'^^ cxp}-^ f dr f d'^r *%y(*2)-c2$. V^*-2J-$ I, (1) 

with inverse temperature j3 = 1/k^T. The 0(A'^)-invariant potential V assumes the 
form 

l^(*2)=^2^2^A$4^... ^ (2) 

where the factor 1/A^ in the second term ensures a well-defined large- iV limit. For 
vn? > 0, we get an 0(iV)-symmetric ground state, ($) ~ 0, whereas m^ < induces 
symmetry breaking. ($) ^ 0. We can formally split the field into modulus x > and 
direction n with n^ — 1 

^ir,T) = xir,T)n{r,T). (3) 

The functional integration measure then becomes 

S)^*=||>SAr_ii|Dxx'^"'2)'^"'n, (4) 

where ||iSjv-^i|| = 27r^^'^/T{N/2) is the surface area of the unit sphere in N dimensions 
and thus we get for the generating functional 

Zp[J]^\\SN-i\\Js''-'nZf,[n,J] (5) 

with the remaining integration over x (note that n • n = by normalization) 

Z}[n, J]^ J^x cxp IJdTJ d^r ({N - 1) In x - ^ [x^n^ + x^ + V{x^)] + . 
lo 

(6) 

In view of the phase-space factor x^~^ in the measure Q, we get, for A^ ::^ 1, a strong 
enhancement of contributions for large x- On the other hand, the terms in the action 
in ^ such as x^"^ yield a strong suppression at large x- Therefore, the x-integral 
will be dominated by the saddle point of the exponent in © given by 

N-1 .. .-2 IdV 
X 2 dx 

In the large N limit, "h? is the sum of many independently fiuctuating terms and thus 
approximately constant (the same is true for n • V n). Therefore, the saddle-point 
approximation yields a constant field x, i-e., a classical (mean) field 

X^ = $2 « xLss = (*') = 0{N) . (8) 

Due to the steep potential for the field x (in the large N limit), its quantum fluctuations 
can be neglected and we can approximate it by a constant c-number. In the above 
formula, we have assumed a suitable ultraviolet regularization resulting in finite 
expectation values (*2) = o{N) and (*2) = 0{N). 

Consequently, to lowest order in l/iV, we obtain the same results as with a bilinear 
effective action 

ZfiJ] = I S'^^ exp J -i y dr y d^r [t/.^ + ml^4>' - cl^cjy ■ V V - 2i • (/.] I , (9) 
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where the renormahzed constants m^g and c^g depend on the shape of the nonUnear 
potential V. Of course, this hne of argument is very similar to the nonlinear cr-model 
which does also approach a linear theory for TV ^ 1. As one may easily imagine, 
the same procedure can be applied in more general cases. For example, one may 
include arbitrary dispersion relations by adding terms like (p ■ 'V cf) and cp ■ V^^ etc. 
Furthermore, the prefactors in front of the different terms (such as c^) may depend 
on (p , similar to the potential V{(p''). 

The general idea is always the same: a sum of many {N ^ 1) independently 
fluctuating quantities on an equal footing such as is approximately a c-number 
whose fluctuations are suppressed by I/^/N (cf. the de Finetti theorem and the law 
of large numbers). Therefore, we shall use the same linear effective action also in 
the time-dependent case (with real time t instead of Euclidean time r) - and with 
a possibly time-dependent (0^)(i), which is due to the growth of the thermal and 
quantum fluctuations (see below). 

3. Phase transition and correlation function 

Having motivated an effectively quadratic action in the large N limit, we can discuss 
the behaviour of the fluctuations during the phase transition, by using a linear 
evolution equation of the form 

f(-V2)0 + G(-V')0 = O, (10) 

where we assumed time-independent coefficients F(— V ) and G'(— V ) initially and 
after the quench. Note, however, that the validity of such linear equations is not 
restricted to large N but a linearization is often also possible in lower dimensions, 
N = 0(1), using different expansion parameters instead of 1/iV. For instance, for a 
spinor Bose Einstein condensate, the gas parameter, measuring the diluteness of the 
system, can be employed, cf. |11) . 

Before the quench, the field is in a stable configuration, i.e., all excitations have 
real frequencies. The field operator can be written in terms of initial modes 

4>{r,t<t,^) = ^|^-^p^A,.£fc,^{e''='^e— '=<*afc,A + e-^'='-e-'=<*aL^}(ll) 

with N orthogonal polarizations, e^-x • £k.\' = ^w- The coefficients Ak depend 
on the functions F{k^) and G{k^) in Eq. PU]) and are chosen such that the mode 
operators a]^ ^ ^^'^ ^k,\ creating or annihilating a quasi-particle of wavenumber k and 
polarization A obey the usual bosonic commutation relation [0,^ x, 0,^1 x'l ~ ^kk'^'xx f^^ 
a bosonic field <p. (In the following, we will take the basis vectors Sa labelling the field 
components as the polarizations £k,\-) 

After quenching the system through a (second-order) phase transition, the 0{N)- 
symmetric state is no longer the ground state of the system. As a result, some modes 
become unstable and thus the (instantaneous) dispersion relation from ^ dips below 
zero for some wavenumbers, w^ ^ < 0, cf. Fig. [31 The time-dependence of cj) can be 
expressed in terms of instantaneous frequencies 
r fjNi. 
0(r, t > tout) = J T^^Tl^e''^" [Bfce.-'"-'* + Cfee»"°-*] ek,xak.x + H.c. (12) 

with initial mode operators ak,\ and complex coefficients Bk and Ck depending on 
the shape of the transition. Note that many of the modes are still stable after the 
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quench and can be viewed as "virtual" short-lived defect-antidefect pairs annihilating 
immediately after their creation. These stable modes form a time-independent 
contribution to the correlation function. The field linearization outlined in the previous 
section (in the large N limit) is still valid for a finite time after the quench - as long 
as the relative contribution of unstable modes is small compared to the mean field 
(02) = OiN). 

Bearing in mind that we start in the 0{N) symmetric phase, which has (0) = 0, 
we cannot simply insert n = (</>)/ 1(0)| into the general expression for the winding 
number given by Eq. pop below. Instead, in order to address the winding number 
statistics, we must employ some appropriately defined direction operator fi, see |30j . 
Since the new phase grows from unstable fluctuations, the basic domain structure 
with broken 0{N) symmetry and thus also the topological defects are seeded by the 
field distribution of the unstable modes shortly after the quench. Stable excitations, 
on the other hand, can be interpreted as virtual defect-antidefect pairs coming only 
briefly into existence. Especially for large k, these short-lived virtual defect-antidefect 
pairs would dominate the defect number in a given volume for any fixed time due to 
the phase space factor k^^^ stemming from the d^fc-integral. This corresponds to 
the strong singularity of the two-point function l/|r — r'\^^^ in A^ ^ 1 dimensions. 
Since we are not interested in these short-lived pairs, but rather in real long-lived 
defects originating from unstable fluctuations, we introduce a time-averaged direction 
operator 



i I dtg{t)cj>ir,t), 



(13) 



with (temporal) regularization function g{t). 

This time average supresses all short-lived pairs at large k. Since the spatial 
behaviour of the n-correlator P^ will be dominated by one wavenumber fc,, see 




n U 



fc* 



Figure 3. Two generic examples for the evolution of the dispersion relation 
following from Eq. I|10|l during a symmetry-breaking phase transition. Initially 
(dotted line), all fc- values are stable, uj'^(k) > 0. At the critical point (dashed 
line), the dispersion relation touches the fc-axis, and after the transition (solid 
line), modes in a finite fc-interval become unstable, ui'^{k) < 0. The left panel 
corresponds to a case, where uj'^{k = 0) = m? in Eq. ((Sjl remains positive while c^ 
changes sign; while in the right panel, u}'^{k = 0) = m^ becomes negative. In both 
cases, however, there is a dominant wave vector fc* (for large N), as indicated by 
the vertical arrows. 
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Eq. (J18p below, and the time-dependence cancels due to the normalization of n, our 
result is independent of the choice of g{t). The only conditions on g{t) are that its 
support lies well within the region of applicability of the effective action ^ and that 
g{t) is smooth enough. More precisely, the spectral width of g{u!) determines the 
uncertainty of A: = fc* ± Afc which should be small Ak <C k* . 

Via the same arguments as in the previous Section, the normalization factor 
Z = 0{N) can be approximated by a c- number for large N 



Z^ 



dtg{t)4>{r,t) 



2 



1 2\ 



dtg{t)^{r,t) 



+ 0{N). (14) 



Obviously, the rapidly oscillating parts of </>, i.e., stable modes, are smoothed out 
in Eq. (|13p and only the exponentially growing excitations remain. Note that the 
leading order of the normalization factor, Z Ri Z , would yield exact normalization of 
the expectation value, (n^) = 1, but only approximate normalization of the operator 
itself, n^ « 1. Proper normalization of the direction operator, -h? = \, can only be 
achieved by including operator- valued higher orders of Z as well. In view of the mode 
expansion (|12p. we get for the time-averaged direction operator 

1 /■ d^k 

^'='-Cfc£fe>fc,A + H.c. , (15) 



Z J (2^)^/2 

unstable 

where the integration domain in k space runs over all unstable modes. The coefficients 
Cfc include the above time-averaging and the exponential growth of the unstable 
modes, as well as the original coefficients Bj, and C^, and possibly the initial 
distribution, e.g., of thermal origin, of the occupation numbers of the different modes. 
The calculation of the correlation function is now straightforward. Taking the 
basis vectors Ba labelling the components (pa of the field as the polarizations Ek,\ and 
approximating the normalization factor by its expectation value, we obtain 

{na{r)n,{r')) = || | ^e^M-^')|(7.p . (16) 

In the case of spatial homogenity and isotropy (as opposed to internal isotropy of 
the N fields 0a) after the quench, where no direction is distinguished, we can adopt 
spherical coordinates and carry out the angular integration 

(".(-)».(r')) ^ ||i^/.«~-.|C.ri<t^, (17) 

where the index v of the Bessel function Ji, is given by i^ = N /2 — 1. The remaining 
k integral is often dominated by one wavenumber fc*: In large dimensions A'' ^ 1, 
the phase-space factor k^^^ strongly favors contributions near the largest unstable 
wavenumber, cf. Fig. [31 whereas for small N = 0(1), the dominant wavenumber 
is usually determined by the largest coefhcient Cfe, i.e., the fastest-growing mode. 
In either case, as a result of the existence of a dominant wavenumber fc.^ and the 
ensuing application of the saddle-point integration, the correlation function of the 
time-averaged direction vector now reads 

r ( \^ ( '^\ A 2'-r(;. + l) J.(fc,|r-r^|) 

{na{r)n,{r )) = 5^, ^ ^k^r - r'\Y ' ^^'^ 

Note again that this is only a part of the direction correlations, generated by the 
unstable modes, and thus responsible for the formation of topological defects. All 
stable modes have been averaged out in definition p^ of the direction operator h. 
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For large iV ^- oo, we can use the well-known asymptotic behavior of Bessel 
functions )32| and obtain Gaussian correlations 

(-a(r)n,(r')> = %i exp ( - MlZLzZ:!! \ , (19) 



N "^ \ 2N 



which fall off on a typical distance y/N/k^,. Note that, since the corrections to the 
normalization factor Z of the direction operator are suppressed for large N, cf. ([M]), 
these Gaussian correlations represent an exact result. The typical oscillatory behavior 
of Bessel functions is exponentially suppressed for large A'' — )■ 00. For finite N, the 
spatial decay of the correlations is too weak such that the oscillations are retained. 
However, for finite N, one should also be aware that corrections to the normalization 
Z would become important, though the leading order of the directional correlations 
would still be given by expression P^. 

4. Topological Defects 

4.1. Winding numbers in dimension N: Classification of topological defects 

As a quantum version of the classical winding number counting the net number of 
topological defects within a domain A^ in A^ dimensions |33| , we introduce a winding 
number operator *Jl for the direction operator n of the 0{N) field, cf. [30] 

P pa/37... r 
MM) = /^'j^-^-\\Sm-i\\ I dS-n%dpn'){d,n-)... , (20) 

dM 

where ||5Ar_i|| ~ 2tt^/^ /r{N/2) is the surface area of the N dimensional unit sphere. 
The operator 9T (and similarly the classical winding number) remains unchanged if we 
subject the fields n to a global internal 0{N) rotation which can be smoothly deformed 
to the identity transformation. Improper rotations, which are also 0{N) symmetry 
transformations (but cannot be smoothly deformed to the identity transformation), 
generally do not leave 9T invariant. 

To exemplify this, let us consider the simplest nontrivial case, N — 2. 
Parametrizing the field direction through a single angle, n^ — cos 6 and n^ ~ sinfl, we 
obtain 

m ^ Y idi-ve, (21) 

c 
where dl is a tangent vector along the integration contour C. A rotation of the field 
direction n' = D{(f>)n by an angle (j> will only add a phase to the angle operator 6 
and the winding number (|2T1) is unaffected, see also Figure |4l Besides rotations, the 
group 0{N) also contains improper rotations, i.e., combinations of a rotation with a 
reflection, which do not conserve winding number. This leads us, apart from the well- 
known monopole or spin-vortex configurations depicted in Fig.[4l to field distributions 
with a nontrivial winding configuration, where the field points away in one direction 
and towards the defect core in the other, see Fig.[Sl 

This can be generalized to higher dimensions. For example, in three dimensions, 
any rotation of n can be parametrized by three Euler angles, i.e., as a sequence of three 
two-dimensional rotations. Obviously, inversion of an odd number of field direction, 
e.g., n — >■ — n, cannot be expressed in that way but corresponds to improper rotations 
again so that the winding number picks up an additional factor —1. Hence, the typical 
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Figure 4. Typical field configurations for topological defects in two dimensions. 
The shown field distributions (hedgehog, left vortex, inverted hedgehog, and right 
vortex) can be transformed into one another by (internal) SO{2) rotations, which 
preserve winding number. Hence, all four belong to the same homotopy class and 
have the same windung number 0^ = +1. 



/tt\ 

i: • r 
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Figure 5. Antidefect with 9^ = —1 ("hyperbolic defect") which follows from any 
of the plots in Figure|4]after field inversion in one direction, i.e., after an improper 
rotation. 



0^ = ±1 configurations liave the field pointing either away from (+1) or towards (—1) 
the defect core (in all directions). Other nontrivial field distributions, e.g., those with 
n pointing away in some and towards the core in all other directions, can be obtained 
through rotations. The strikingly different appearance of positive and negative defects 
in two and three dimensions is typical for even and odd number of dimensions A^: Field 
inversion, n -^ —n, will leave the winding number *Tt -^ (—1)^91 invariant in even 
dimensions while for TV odd, VI picks up a factor (—1)^ = —1. As we will see in 
an example below, this has implications for the net defect number created during a 
quench to the broken-symmetry phase. 

4-2. Winding number variance 

Since, for a fully 0{N) invariant initial state, the probabilities for creating positive 
and negative defects are equal, the expectation value (DT) of the winding number 
operator (j20p is zero, but the variance (0^^) is generally non- vanishing and can be used 
to quantify the net number of defects created during the quench. Due to the totally 
antisymmetric tensor Sabc... on the right hand side of Eq. (j20p . each field component fia 
must appear exactly once in *Tt. Hence, if different field components are uncorrelated, 
cf. Eq. dH]) 

{n''{r)n\r')) = f{r - r')5''' = /(|r - r'|)<5"^ , (22) 

the expectation value (W) must break down into the product of iV two-point functions. 
The first identity in ([^^ holds when the system is homogenous, while the second is 
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valid when the system is isotropic as well. This particular form of the two-point 
function {n"' {r)n}' (r)) can be inferred from symmetries: Although the phase with 
broken 0{N) symmetry is neither homogeneous nor isotropic, the expectation values 
still are, because we start from an homogeneous and isotropic phase and any realization 
of the symmetry breaking is equally likely. 

With (|22p , there results the following expression for the winding number variance 



{^'{M)) = „^^ ,„ </ dSa.e'^^-'- <f dS'^,e^' "'-''■■■ 
WSn-iV J J 

dM dM 

X [f{dpd'p,f) -{N- l){dpf){d'p,f)] {d,d'^,f){dsd's,f)... , (23) 

where we exploited the symmetry under exchange of index pairs (a, a) O (&,/3), and 
also changed the order of taking spatial derivatives and expectation values. Due to the 
total antisymmetry of the integrand, this expression can be regarded as a collection of 
surface integrals over {N — l)-forms (in the unprimed indices while keeping the primed 
indices fixed and vice versa). It becomes therefore possible to apply the generalized 
Stokes theorem, f^ dtj — Jgj^ w, and rewrite (|23p as a volume integral over iV- forms 

i^HM)} = -^ ^e"^"-e" ^ ^ - / d^rd''r'{d^d'^,f){df,d'^,f){d,d'yf)... . (24) 

M 

This expression is generally valid provided different field direction are uncorrelated, 
{hafib) oc 5ab, cf. ([22|) . The next step is less obvious, but again we can take advantage 
of the inherent symmetries which enabled us to write the two-point function (|22|) 
through a function depending only on the separation L = |r — r'|. Using the identity 

do^d'^.f = {dcd'^^L'-)^ + (a„L2)(a;,i2)|V (25) 

the derivatives under the integral can be evaluated. In view of the total antisymmetry 
of the Levi-Civita symbols, all terms containing more than two first derivatives of L^, 
e.g.; {daL'^){d'^,Lp'){dpL^){d'g,L^) cancel after summation. It follows 

(^^) = p^^"''^-e"'^'-'-/d~rdv|(a„a;,L2)^ + iv(a„L^)(a;,,i^)0 

N-l 



x{d,d'^,L^){d,d'yL^)...(^^^ . (26) 



So far, we did not specify any particular coordinates, but in order to simplify this 
expression, it is advantageous to adopt Cartesian coordinates for which dad'^,L^ ~ 
—2Saa'- After performing the summations over a, a' etc. and with df/dL^ = 
{l/2L)df /dL, we finally obtain 

M 

which holds for arbitrary integration manifolds M. and any dimension N provided 
the two-point correlator has the form (|22|) . i.e., the correlator {h'^{r)h^{r')) depends 
only on the separation between r and r'. If, however, homogenity or isotropy are 
broken but different directions are still uncorrelated, the most general expression for 
the winding number variance is given by Eq. (|24p. 
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4-3. Spherical volumina 

The 2N integrals appearing in ([?7|) are, for arbitrary volumina, generally difficult 
to evaluate directly, but can be simplified drastically by exploiting symmetries of 
the integration manifold Ai. In particular, we will demonstrate that for a sphere 
of radius R - an isotropic manifold defined through Ai = {r : r^ < R^} - Eq. 
(|27p reduces to a single integration over one remaining variable. To this end, let us 
instead of the winding number variance (tn^) consider its change with radius R of the 
integration volume. Taking the derivative of (j27p with respect to R and using spherical 
coordinates, one of the radial integrations will break down to yielding a factor R^^^ 
from the integral measure [note that the integrand of (j27p is symmetric in r and 
r']. One of the angular integrations d^~^n' can be performed due to isotropy of the 
manifold Ai and gives the surface area |j57v-i|| of the N dimensional unit sphere. For 

we get 

NN\ „„_, f . f .„_,„ Jn d 



the radial derivative of (*Tl^) 



OR 



l^AT-l 



tR 



N-1 



dr 



iN~l 



Q. 



L^-i dL 



Of 
dL 



N 



(28) 



The remaining volume integral d r = Jn drd H} with the Jacobian J^ can be 
tackled by a coordinate transformation r -^ L — r — r' . Adopting spherical 
coordinates for Z,, we see that the (radial) measure L^^^ just cancels the factor 
\/L^~^ such that the remaining L dependence is just a partial derivative and the 
modulus integral dL can be performed easily. Regarding the angular part, we exploit 
isotropy of A^ once more and choose the angles 0i...0Ar_i such that the separation 
L = 2i?sin0i depends only on 0i. After substituting y ~ sin0i and performing the 
integrations over the remaining angles, we then further get 

1 



dR 



= 2NNl 



\S. 



N-2\ 



\\Sn-i\ 



R 



N-l 



Jdy^l^ 



-N-3 



V 



-i<-') 



N 



.(29) 



tN~3 



-N-1 



Using the identity ^1 - y^ = -{l/N - l)y^<9(\/l - 2/^ /y^-^)/dy, Eq. ^ 

can be rewritten and, after integration by parts, the change of the winding number 
variance reads 



OR 



NN\ \\Sn- 



N-2\ 



N-1\\S. 



N-l\ 



dy 



V^ 



-N-1 



y 



1 d 



V 



N-l 



Rdy 



(Ry) 



N 



dR 



my) 



N 



(30) 



In the factors right of the derivative, the integration variable y appears only in the 
combination Ry. Hence, using the chain rule, d/dy can be replaced by {R/y)d/dR 

and it follows 

1 



9(91^) _ iViV! ||5jv_2|| d 
OR ^ iV-lT^ 



N~l 



dR 



R" JdyV^ 



-N-1 



y" 



-li^'^y^ 



N 



(31) 



Both sides of this equation are partial derivatives with respect to the manifold radius R 
and can be formally integrated, where the constant of integration follows by requiring 
{'Vl^{R = 0)) = 0. Finally with y = sin(6'), we can simplify the result to the more 
compact form 

tt/2 

m 



(fn^(i?)) = —R 

IT 



N 



de 



— cos 



^|(2i.sin< 



N 



(32) 
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i.e., the winding number variance can be determined through a single angular 
integration. For short-ranged correlations, where df /dL falls off fast enough at large 
distances L 

'^<OiL'^):^>^, (33) 



a surface scaling |1H 130] 

(m'iR)) oc i?^-i (34) 

follows for large R: For sufficiently large R, the main contribution to ([5^ arises from 
small 9. By formally extending the integration to infinity and approximating cos w 1 
and sin 6^6, one can see that the integral is proportional to 1/R and we obtain scaling 
of the winding number variance with the surface area of the enclosed volume. Note 
that expression (|32|) and thus also the surface scaling for short-ranged correlations 
(|33p hold in any dimension N as long as the direction correlator takes the form (22) . 
In particular, the expression for {'^l'^) derived in |11) will be reproduced for N — 2. 

We stress that our calculation based on very general principles predicts a 
surface scaling law for the net defect number variance. This surface law should 
be contrasted with the volume scaling law obtained from another frequently used 
and phenomenologically motivated approach |2]. If one would randomly distribute 
hedgehogs and antihedgehogs in a given TV-dimensional space with the same fixed 
densities, the total number of defects (hedgehogs plus antihedgehogs) inside a 
given volume would scale with the volume itself. Their difference (i.e., hedgehogs 
minus antihedgehogs), however, would scale with the square root of that volume 
(corresponding to a one-dimensional random walk of the winding number) and thus 
the net defect number variance scales with the volume. 

As far as we are aware, the rather generic result obtained here, together with our 
previous results in [TTl [5D], represent the first proof by direct calculation that (the 
quantum version of) the Kibble-Zurek mechanism leads to a surface law. 

4-4- General convex volumina 

The surface scaling for short-ranged interactions can be interpreted as a variety of 
a "confined" phase, where bound defect-antidefect pairs, typically separated by a 
correlation length, occur. Pairs entirely within or outside the manifold Ai will not 
contribute to the net winding number and thus also not to (9T^). Only those pairs, 
where one of the partners is inside and the other outside will contribute. Since it is not 
clear whether the defect or the antidefect is contained within Ai, calculation of the 
winding number would be similar to a random walk with number of steps proportional 
to surface area. From this simple picture, it can be anticipated that the surface scaling 
is not restricted to spherical manifolds but might apply in any situation where the 
surface of the manifold is sufficiently "smooth". 

In the following, we will thus consider arbitrary convex manifolds and show that 
the surface area scaling of the winding number variance is retained. We parametrize 
the convex manifold in spherical coordinates 

M = {r:r^< R^a^{i}i, ...^n-i) = R^a'^i^)} 

with sri + {1 ~ s)r2 e V V ri, ra € T^, s e [0, 1] (35) 

through a scale factor R and a shape function a(ri) = 0(1), similar to a sphere, where 
we would have a{Cl) = 1. Obviously, the surface area DM oc R^~^ so that we need to 
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show (9T^) oc R^~^ as well. For convex manifolds, the volume element in expression 
(|27p for the winding number variance reads in spherical coordinates 

_Ra(f2) 

^r^ fd^-^njN{n) I drr^-\ (36) 



where the Jacobian JN{n) is taken at r = 1. As for spherical volumina, cf. ((28|). the 
scale factor R appears only in the upper bound of the radial integration such that we 
can use the same trick as before and calculate d{V(^) /dR instead of the variance itself 

However, in contrast to the spherical case, we cannot simply perform the angular 
d^~^il integration because of broken isotropy (of the manifold): The separation 
L = \r — r'\ generally depends on both azimuthal positions ft as well as Jl' and 
not only on their relative angles. The coordinate transformation r ^ L — r r' and 
subsequent integration in spherical coordinates yields 

^^ =2^^i?^-iy'd^-il7'dOJ^(f2')^^(f^)a^("') (-|^ [L„ax(n,f^')]) , 

(38) 

where Lmax is the distance from a point $7' on the surface dM. to the opposite 
boundary of M in direction Q, and is determined through L^ax = |'"max[^(i^)] — 
rj^jj^(ri')|. [At this point, it should be mentioned that one needs to be careful regarding 
the azimuthal positions: $7' and $7 are those of r' and r as seen from the predefined 
centre of the manifold, cf. ([35]) , whereas fJ is the azimuthal position of r as seen from 
r'.] A scaling of the winding number variance with surface area is obtained if the 
integral in (|38|) is proportional to 1/R for sufficiently large R - which can indeed be 
expected for short-ranged correlations ([55)) because r„iax and r'^^^^ both scale with R 
and thus L^ax = R^ma^{^, ^') « R. 

We evaluate the scaling behaviour of ([55)) by choosing suitable coordinates for 
the integration: In view of the convexity of M, the angles Ct, can be chosen such 
that for each Jl' there exists at least one angle 0i of Ct for which the separation 
imax = -RAmax chaugcs from zero to order R while keeping all other angles 02 to 
4>N~i fixed at arbitrary values. The 0i integral appearing in ()38p will be performed 
first, it reads 



■ JV-2 

sm ( 



df 1^ 

-^[i?A„,ax(</'l)] . 



(39) 



Obviously, the integrand vanishes for all (pi with A,jiax('/'i) = 0. It becomes thus 
possible to change the domain of integration to [(/)_, 0+] C (0,7r) where A„iax(</') > 
for all 4> G ((/>_, 0+) and A,„ax('/') ~ OV0 ^ ((/>_, 0+). From convexity of the 
integration manifold M also follows that the maximal separation Amax cannot have 
(local) minima in (0_, </>+), i.e, Amax nrust be monotonically growing in the vicinity 
of 4>- and decreasing near (/>+. Furthermore, by demanding short-ranged correlations, 
(|33p. df/dL becomes arbitrarily small for large separations. (See, by contrast, the 
case of long-range correlations we discuss in the following section). This means that. 
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for sufficiently large R, the main contributions to integral ([5^ stems from regions 
where A,„ax is small but non-zero, i.e., near 0_ and 0^_. and we can approximate pop 



dcpi sin 



Af-2 



N 0a+« ^ 



with small e ^ 1. For all angles (p E (</)_ + e, 0+ — e), the distance A^ax > but 
the integrand is negligibly small due to short-range correlations ([55)1 . In the integrals 
on the right hand side of PO)) . we can assume that Amax is strictly monotonic and 
analytic such that Aniax('/'i) is bijective in (0_ , ^_ -I- e) and similarly in {(j)^ — e, </>+). 
It becomes therefore possible to change the integration variable from 0i to A,„ax and 
we finally obtain for the 01 integral ([39| 



7r 



6i sin^- 



df 

-^[i?A„,ax(</'l)] 



TV 



(41) 



a-± Q 

where we have extended the integration to infinity because, due to short-range 
correlations, the integrand becomes negligibly small for sufficiently large R. In this 
limit, the integral yields a constant (independent of system size R but of course still 
depending on all other angles (j)2...(j)N-i as well as J7') and the scale factor R appears 
only as prefactor 1/R. We are thus led to the scaling of ([55)) with R^^^ for large R 
and thus, after integration, obtain for the winding number variance 

(<yi2(i?)) oc i?^-\ (42) 

i.e., the desired surface scaling result. It should be noted that the R^~^ dependence of 
d{^) /dR derived in this section is only valid for sufficiently large R. When integrating 
d{^) /dR to obtain the scaling behavior of the winding number variance, regions with 
small i?, where d{W)/dR might behave differently, are also included, and potentially 
contribute subleading orders to the winding number variance, which scales with surface 
area for large R. 

5. Deviations from the Surface Law 

We have obtained the surface scaling ([5^ for spherical volumes in high dimensions, 
and generalized this result to arbitrary convex volumes ([^^ . The salient ingredients 
of the derivation have been the short-range nature of correlations, Eq. ()19p . and the 
symmetries of the correlation function ([22)) . By relaxing either or both of these 
conditions, and depending on the precise functional form of the two-point correlator, 
it is then possible, as we will now show, to obtain various interesting scaling behaviors 
with system size. 

Let us consider broken spatial isotropy of the correlation function ([^^ . e.g., 
through an externally imposed lattice. (Note that {nafib) oc Sab is still isotropic in 
the field directions.) Explicit examples for periodic and short-ranged correlations will 
be discussed in subsections 15.11 and 15.21 To remain general, we assume that in ([22)) 
only the first equality holds, and that the (real space part of) the correlator separates 
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into the product 

N 

f{r-r')^l[g^{r^~r'J, (43) 

a=l 

where the correlations g,^ in different directions can in general have a different 
functional behavior. We insert this expression into (j24p . and consider the volume 
enclosed by a A^-dimensional box with size Li x L2 x L3.... The second-order partial 
derivative of the correlator (|^5)l gives 



9'a9L^ 


1 5aa' 




9a gay 





9a 9a' I ( 9a Ha Ha \ x ('44'1 



dad'^.f = -f 

_9a9a' \9a 

where a prime denotes the derivative, e.g., g^ — dga{x)/dx. Inserting 
the above expression into ((24|). we see that any terms involving products 
{g'al 9a){g'ai / ga'){g'j3/ gi3){g'i3, / gp') etc. cancel during summation due to antisymmetry 
of the e pseudotensors and we have 

I \\Sn-\V J J ' \9a 9a 9a 

J?l_?k?k)(?l^?ki).,^, (45) 

\gp 913 913 J \9-i g^y g-(J 

This correlation function can be evaluated further by accounting for the broken 
isotropy of / and regarding the integral over each each spatial direction independently. 
It follows 

Li3 Lj3 

\N-1atat\ r r 

JV-l // N~2 I I 

/35^a 



^'> = ^^lljj2''' Y.y9^A^)-g^aiLa)] n jdr.jdr', 



9/3 9 13 ~ 9 13 gpgp 



(46) 

So far, this result is general. We now are going to illustrate it with a few concrete 
examples for the correlation function. 

5.1. Periodic Order in the Direction Operator Correlations 

As an example for a nonlocalized correlation function, we consider a purely periodic 
functional behavior, which can, for example, be derived within the following simple 
model. Using the ansatz for the fluctuating order-parameter field 

<Pa{r) = ]^cos(A:a::6 -f ipab) (47) 

b 

and averaging over N x N independent fluctuating angles (pab S [0,27r], one obtains 

2tt 

vanishing expectation value {(f>a{r)) ~ Ylb 'Tn I d'Pa.b cos{kxb + fab) = 0. The two- 


point correlations are then calculated to be 

{cj)a{r)Mr)) = l6abY[cos{k[xc~x',]), (48) 

c 

and explicitly break the isotropy of the system by imposing periodicity along all N 
spatial axes. 

The integrals in (|46p now grow with system size La for an even number of 
dimensions (the integrand is strictly positive for even N), whereas the first term in the 
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sum, 5^(0) — da i-^a) ~ 1 ^ cos^(fcLfj), is a periodic function in the system size L^- 
Thus the correlation function (|^5|) implies a scaling of the winding number variance (for 
even N) with the square of the surface area of integration volume F = ii x _L2 x -^3--- 

mA ^ y^[i- cos'' (kL^)]^. (49) 



even 



The factors 1 — cos^ (kLa), however, will still lead to vanishing winding number 
variance if the size La = m:/k matches the periodicity of the correlations in all spatial 
directions. For an odd number of dimensions, on the other hand, the integrands of (|46p 
oscillate periodically around zero and no scaling of (9^^) with system size is obtained, 

(fn2) = oiv°). 



\\ o \ o \ o \ 

t r ^ o t o I 

I c t o i o t 

t c i o t o i 



Figure 6. Distribution of topological defects in the presence of periodic order in 
two dimensions. The unit cells of the lattice are depicted as gray and white areas, 
while the defect cores are depicted as yellow for 91 = +1 and green for ?U = — 1, 
cf. Figs. [4] and [5l 

This entirely different scaling - surface squared with oscillating prefactor for even 
dimensions versus essentially no net defects in odd dimensions - follows directly from 
the behavior of the winding number under field inversion described in Sec. 14.11 (see 
also Figs.|4]and[5]): Due to the strict periodicity of the correlation function (|48| . the 
field orientation must be inverted after half a unit cell n{r) = —n(r + eair/k) in every 
direction Ca, see Fig.[6l Hence, for each defect with topological charge Vt = +1, there 
sits another defect with charge VI = (—1)^, a distance ir/k apart in every lattice 
direction. This means for even dimensions, one has alternating layers of defects and 
antidefects, and thus the surface squared scaling with oscillating prefactor results, 
whereas in odd dimensions, each layer contains positive and negative defects, so that 
the total winding number essentially vanishes for a box-like integration manifold. 
(However, if the box-shaped integration manifold were tilted with respect to the lattice, 
a non- vanishing net defect number could occur in odd dimensions as well.) 

5.2. Localized (Gaussian) Order 

As a second example, we take short-ranged correlations, e.g., a Gaussian correlator 
similar to (fT9|l derived for spatial isotropy in large N. Now, each of the double integrals 
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in (j46p scales with i^, such that the winding number variance is proportional to the 
surface area of the enclosed volume 

^V ^ n Lp h^ [g^p~\L)g';,{L)-g^p~^g'p{L)g'p{L)\ . (50) 

The above relation implies that the surface scaling found for spatially isotropic short- 
range correlations, cf. Eq. ([32]), is shown to be preserved also for the case of broken 
spatial isotropy. This should be contrasted with the result of the previous subsection, 
where it was shown that for long-range periodic order, no scaling proportional to the 
surface area is obtained. 

6. Conclusion 

Starting from a rather general 0{N) symmetric action, we derived a quadratic 
effective action for the field fluctuations, and studied the instability of modes during a 
phase transition from the initial 0{N) symmetric ("paramagnetic") phase to a broken 
symmetry phase in which "ferromagnetic" order develops. Since we are interested 
only in the exponentially growing modes, which dominate the development of the 
domains and thus also of topological defects, we introduced a systematic averaging 
procedure for the fluctuating field operator directions. We have demonstrated that 
the direction field correlation function can be explicitly evaluated. Our results are 
very general: The initial state (such as temperature, dispersion relation) enters the 
coefficients Bk and Ck in (fT2|) only. Thus, while the initial state affects the field 
correlator {(l)a{t,r)(pi,{t,r')), it leaves the direction n-correlator P^ invariant. The 
external conditions after the transition determine the final dispersion relation which 
yields the value of fc*. Apart from a simple rescaling of this value, our main results 
for the homogeneous and isotropic case (surface scaling law etc.) thus do not depend 
on the initial and the final state. 

This correlator of the field direction then directly determines the statistics 
of the net defect number created by the quench. A general expression for the 
winding number variance was derived and was shown to scale with the surface area 
of the enclosed convex volume in any (large) dimension N. Such a surface area 
scaling follows from general arguments based on the assumption of short-ranged 
correlations together with the internal 0{N) and spatial symmetries (homogenity and 
isotropy). For large enough samples, it should be possible to verify this prediction 
experimentally. The detection method depends on the specific experimental setup: In 
Bose-Einstein condensates, the topological defects are phase or spin vortices (in two 
spatial dimensions) which can be imaged either by absorption (resolving the core as a 
hole in the condensate |25)). or using polarization-dependent phase-contrast imaging 
|34| . In superfluid 'fte, vortices can be detected with NMR techniques [3^1 • Finally, 
in conventional condensed matter systems with magnetic order like ferromagnets, 
topological defects can be detected for example with Lorentz transmission electron 
microscopy [36| . 

If the correlations are not short-ranged, the scaling law can be very different: 
E.g., for purely periodic (i.e., long-ranged) correlations on a lattice, the scaling law 
strongly depends on the parity, i.e., whether the number of dimensions N is even or 
odd. In even dimensions, the winding number variance basically scales with surface 
area squared with oscillating prefactor. By contrast, for odd N, the net defect number 
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inside a box of size V — Li x L2 x L3... effectively vanishes. This entirely different 
behavior can be explained by the strict periodicity of the field directions together with 
properties of the winding number operator, i.e., how positive and negative defects 
are aligned on a lattice: When inverting the field direction, n — > — n, the winding 
number 91 acquires a prefactor (—1)^. Hence, configurations where the field points 
towards to and those where it points away from the defect core must lie in the same 
homotopy class {^ = +1) for even dimensions and in different classes (91 = ±1) in 
odd dimensions. The strict periodicity of the correlation function, on the other hand, 
implies a field inversion after half a period (along any of the lattice directions). Thus, 
one has alternating layers of positive and negative defects for even dimensions, while 
for odd dimensions positive and negative defects appear in each layer. 

Further possible directions of research include the extension of our results for the 
winding number variance to finite, experimentally accessible values of A'^. We stress in 
this context that the replacement of the quantum normalization factor Z introduced 
through the time-averaging of the direction vector, by its classical expectation value 
in the correlator (|16p . remains the only finite N approximation made in our analysis. 
This statement holds under the condition that we have a system for which the effective 
action for the field remains still quadratic for finite N (like in a spinor Bose-Einstein 
condensate) leading to a dispersion relation with a dominant wavenumber fc* > 0. The 
finite N extension can then be achieved by evaluating Eq. ([T5)) to the next order(s) in 
1/N^ i.e., by calculating the corrections to the operator- valued normalization factor 
(|14p entering the average quantum direction vector in ([TS]) . 
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